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1  SUMMARY 

Analysis  of  time-series  data  can  yield  remarkably  accurate  estimates  of  the  frequency  of  a  satellite’s  brightness 
modulations.  These  “apparent”  or  synodic  frequencies  can  vary  in  time,  differing  from  the  actual  rotation  rate  of  the 
object  by  an  amount  that  depends  on  the  relative  angular  motion  between  the  satellite,  illuminator,  and  observer  for 
reflected  light  measurements  (or  between  the  satellite  and  observer  for  thermal  emission  measurements).  When 
detected  with  sufficient  accuracy,  such  synodic  frequency  variations  can  be  exploited  to  characterize  an  object’s 
rotation  state,  using  an  analysis  that  does  not  require  any  a  priori  knowledge  of  the  object’s  shape.  For  instance,  this 
shape-independent  analysis  method  can  be  used  to  derive  spin  axis  orientations  and  sidereal  rotation  rates  for 
spinning  objects.  Remotely  determining  such  rotation  parameters  can  be  useful  in  many  circumstances,  such  as 
when  performing  anomaly  resolution  for  satellites  that  have  lost  stabilization.  Unfortunately,  synodic  variations 
cannot  be  detected  by  ground-based  observers  for  many  objects  due  to  low  rates  of  relative  angular  motion.  This  is 
especially  true  for  non-specular  objects  in  deep-space  and  geosynchronous  orbits.  In  these  cases,  deriving  spin  axis 
orientations  can  be  accomplished  using  a  shape-dependent  method  that  employs  a  model  of  the  shape  and 
reflectance  characteristics  of  the  object.  Our  analysis  indicates  that  a  simple  cylinder  shape  model  can  suffice  to 
characterize  rotation  states  for  upper-stage  rocket  bodies,  although  even  this  relatively  simple  model  requires 
significantly  more  computation  than  the  shape-independent  approach. 

2  INTRODUCTION 

Ground-based  optical  and  radar  sites  routinely  acquire  resolved  images  of  satellites,  yielding  a  great  deal  of 
knowledge  about  orbiting  spacecraft.  However,  the  important  population  of  deep-space  satellites  often  cannot  be 
resolved,  because  they  are  too  far  away  or  too  small.  These  must  be  characterized  using  methods  other  than 
imagery.  In  this  regard,  temporal  photometry  (i.e.,  measurements  of  whole-object  brightness  as  a  function  of  time) 
can  be  very  valuable.  For  instance,  satellite  rotational  motion  can  be  detected  and  characterized  by  analyzing 
periodic  brightness  variations  [1-5].  Temporal  brightness  patterns  can  also  be  exploited  to  characterize  stabilized 
GEO  satellites  [5-10],  The  objective  of  this  paper  is  to  outline  how  temporal  photometry  can  be  used  to  characterize 
non-resolvable  object  rotation  states,  focusing  on  objects  well  above  low-Earth  orbit  altitudes. 

3  OBSERVATIONS  OF  NON-RESOL VED  ROTATING  OBJECTS 

This  analysis  focuses  on  broad-band  temporal  photometry  of  Earth-orbiting  satellites  (i.e.,  time-series  measurements 
of  whole-object  brightnesses)  measured  by  silicon-based  CCD  detectors  sensitive  to  wavelengths  in  the  0.4  to  1.1 
pm  range.  In  this  spectral  range,  reflected  sunlight  tends  to  dominate  satellite  brightnesses,  especially  for  deep- 
space  objects  where  contributions  from  earthshine  and  other  sources  can  be  neglected  to  first  order.  Upper-stages 
from  rockets  that  have  launched  satellites  into  orbit  represent  one  of  the  most  common  types  of  rotating  object  in 
Earth  orbit.  Extensive  observations  of  such  expended  rocket-bodies  have  been  conducted.  This  analysis  focuses  on 
analyzing  ground-based  observations  of  one  type  of  upper-stage  rocket-body:  Boeing’s  Delta-IV  Stage-2  (D4S2). 

3.1  Delta-IV  Stage  2  Upper  Stage  Rocket-Bodies 

Figure  1  shows  a  rendering  of  a  D4S2  [11],  which  like  many  other  rocket-bodies,  has  a  roughly  cylindrical  shape. 
D4S2s  measure  roughly  12.2  m  in  length  and  4  m  in  diameter,  and  have  been  used  to  launch  a  variety  of 
geosynchronous/geostationary  (GEO)  satellites  as  well  as  several  GPS  satellites  [11,  12].  As  of  2014-07-03,  the  US 
Air  Force’s  Space-Track  web-site  tabulated  orbital  elements  for  eight  D4S2s:  four  for  GEO  and  four  for  GPS 
launches  [12],  The  four  used  for  GEO  launches  now  reside  in  geostationary  transfer  orbits  (GTO),  with 
eccentricities  e  >  0.5  and  mean  motions  of  ~  2  revs/day.  Those  used  for  GPS  launches  reside  in  nearly  circular 
orbits  (e  <  0.03)  but  have  similar  mean  motions  [12], 

3.2  Raven  Telescope  Photometry 

Dining  the  last  few  months  of  2013,  photometric  observations  of  several  D4S2s  were  conducted  by  AFRL’s  small 
Raven  telescopes  located  at  the  Remote  Maui  Experiment  (RME)  site  [13,  14],  The  top  panels  of  Figures  2-4  show 
representative  plots  of  three  of  the  measured  light-curves.  All  of  the  measurements  show  evidence  of  periodic 
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brightness  variations.  The  light-curves  vary  in  duration  from  a  few  minutes  up  to  an  hour  or  longer.  The  longer 
light-curves  often  show  evidence  of  significant  changes  in  the  amplitude  of  the  brightness  variations  (see  Figure  4). 


Figure  1.  Rendering  of  the  Boeing  Delta-IV  Stage-2  upper  stage  rocket-body  (from  [11]). 

3.3  Periodicity  and  Frequency  Spectrum  Analysis 

The  Raven  light-curves  can  be  analyzed  quantitatively  for  periodicity  using  the  frequency  spectrum  and  Fourier 
analysis  methods  outlined  in  reference  [2],  Figure  5  shows  a  representative  periodicity  analysis,  in  this  case  for  the 
first  ~48  minutes  of  the  SCN  29156  light-curve  plotted  in  Figure  3.  The  periodicity  analysis  is  a  multi-step  process. 
First,  a  Lomb-Scargle  periodogram  frequency  spectmm  [15,  16]  is  calculated  (Figure  5,  middle  panel),  and  the 
peaks  in  the  spectrum  above  the  99%  statistical-significance  confidence  level  identified.  Then  a  peak-by-peak 
analysis  is  performed.  This  analysis  subjects  the  original  light-curve  data  to  a  polynomial+Fourier  series  (PFS) 
fitting  process  [see  2],  in  which  each  peak’s  best-fit  period  is  refined  and  then  graded  according  the  quality  of  fit 
provided.  High-quality  PFS  fits  are  identified  with  green  symbols  on  the  frequency  spectrum,  medium-quality  fits 
in  yellow,  and  low-quality  in  red;  the  best-fit  solution  is  identified  using  a  triangle  and  lesser-quality  fits  using 
diamonds  (Figure  5,  middle  panel).  The  final  phase  “folds”  the  data,  or  plots  the  de-trended  brightnesses  against  the 
cyclical  phase  (Figure  5,  bottom  panel).  De-trending  entails  subtracting  the  polynomial  component  of  the  best-fit 
PFS  model.  Cyclical  phases  vary  from  zero  to  one,  and  are  defined  as  the  time  modulo  the  best-fit  period  divided  by 
the  best-fit  period.  Figure  5  shows  two  complete  cycles  of  the  de-trended  data,  following  the  method  commonly 
employed  by  asteroid  astronomers. 
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Figure  2.  RME  Raven  telescope  light-curve  measurement  for  D4S2  rocket-body  SCN  38834  (top  panel),  and 
the  periods  estimated  from  the  data  (bottom  panel)  using  the  sliding -window  periodicity  algorithm  from  [2], 
which  show  no  statistically  significant  periodicity  variations  as  a  function  of  time.  In  this  case,  performing  a 
synodic/sidereal  spin-state  analysis  would  be  fruitless. 

3.4  Sliding- Window  Temporal  Analysis  to  Find  Periodicity  Variations 

The  PFS  fitting  method  also  provides  a  means  to  detect  and  measure  periodicity  variations  occurring  within 
individual  light-curves  using  a  “sliding-window”  technique,  in  which  multiple  analyses  are  performed  for  a  series  of 
overlapping  windows  in  time  sliding  across  the  light-curve  (see  [2]  for  more  details).  Figures  2-4  show  the  results 
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of  sliding-window  analysis  for  the  D4S2  data,  with  the  input  light-curve  data  plotted  in  the  top-most  panel,  and 
estimated  periods  plotted  on  the  same  time  scale  below.  Figure  2  shows  an  example  where  no  statistically  significant 
periodicity  variations  are  detected  (i.e.,  all  of  the  estimated  sliding-window  samples  are  consistent  with  a  flat  line). 
Figure  3  shows  an  example  where  weak  variations  are  detected,  in  this  case  a  modest  but  detectable  downward  trend 
in  the  periods.  Figure  4  shows  an  example  where  stronger  variations  are  detected.  Specifically,  in  this  long -duration 
light-curve  for  SCN  36586,  statistically-significant  period  variations  occur  between  about  1 1:00  UT  and  12:00  UT. 
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Figure  4.  RME  Raven  telescope  light-curve  measurement  for  D4S2  rocket-body  SCN  36586  (top  panel),  the 
periods  estimated  from  the  data  (middle  panel)  using  the  sliding- window  periodicity  algorithm  from 
reference  [2],  which  show  significant  periodicity  variations  near  -  11:30  UT. 
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Figure  5.  Periodicity  analysis  for  SCN  29156,  using  the  first  third  of  the  Raven  data  plotted  in  the  top  panel  of 
Figure  3.  The  top  panel  shows  the  measured  light-curve  in  range-normalized  stellar  magnitudes  (black)  and 
the  best-fit  Fourier+polynomial  series  model  (blue).  The  colored  symbols  in  the  middle  panel  show  the 
statistically-significant  peaks  in  the  spectrum,  and  the  green  triangle  marks  the  best-fit  period  of  341.927  ± 
0.049  s  The  middle  panel  shows  the  raw  periodogram  frequency  spectrum  (blue)  and  the  harmonically- 
extended  spectrum  (black)  [see  2].  The  bottom  panel  shows  the  de-trended  brightnesses  (i.e.,  with  the  best-fit 
polynomial  model  component  removed)  on  the  vertical  axis,  plotted  against  the  cyclical  phase  (i.e.,  the  time 
modulo  the  best-fit  period)  on  the  horizontal  axis,  again  with  the  measurements  in  black  and  the  best -fit 
Fourier  series  model  in  blue. 

4  SHAPE-INDEPENDENT  ROTATION-STATE  ANALYSIS 

Asteroid  astronomers  have  devoted  a  great  deal  of  effort  to  develop  methods  to  characterize  rotation  states  using 
temporal  photometry  (see  [17-19]  and  references  therein)  in  the  absence  of  a  priori  knowledge  of  an  object’s  shape. 
Some  of  these  methods  are  truly  shape  independent,  in  that  they  do  not  require  knowledge  or  inference  of  any  object 
shape  or  reflectance  parameters,  while  others  require  that  object  shape  parameters  be  inferred  from  the  data  along 
with  the  rotation  state  parameters. 

4.1  Single-Axis  Spin  States  vs.  Multi-Axis  Rotation  States 

The  first  step  in  a  rotation  state  analysis  process  entails  discriminating  complex  or  multi-axis  rotation  (precession  or 
precession+nutation)  from  single-axis  rotation  (e.g.,  principal-axis  spins).  Axisymmetric  bodies  in  multi-axis 
rotation  states  produce  brightness  variations  that  generally  require  two  frequencies  to  explain:  one  for  the  spinning 
motion  of  the  body  about  its  principal  axis,  and  another  for  the  precessional  “coning”  motion,  f,/,  [18,  19]. 
Generally,  axisymmetric  bodies  in  long-axis  mode  precession  states  (like  many  asteroids  and  most  rocket -bodies) 
produce  periodogram  spectra  with  the  strongest  peak  at  a  frequency  of  2  (f¥  +  f</i).  Furthermore,  fitting  multi-axis 
rotator  light-curve  brightness  modulations  generally  requires  a  two-dimensional  Fourier  series  [18,  19],  rather  than 
the  one -dimensional  Fourier  series  that  is  sufficient  for  single-axis  spinners  [2].  Fitting  2D  Fourier  series  models  to 
light  curves  can  therefore  provide  shape-independent  means  of  gaining  insight  into  multi-axis  rotation  states. 

Single-axis  spinners,  on  the  other  hand,  produce  light-curves  that  can  be  fit  to  a  high  level  of  statistical  confidence 
using  a  one-dimensional  Fourier  series.  In  this  regard,  finding  that  a  best-fit  ID  Fourier  series  (either  alone  or 
supplemented  with  a  polynomial  series)  fails  to  sufficiently  describe  observed  light -curve  variations  can  be  used  as 
means  to  infer  the  existence  of  multi-axis  rotation,  without  further  characterizing  the  details  of  the  rotation  state. 
This  eliminates  the  need  to  calculate  a  more  computationally-intensive  2D  Fourier  series  fit,  or  the  much  more 
intensive  analysis  to  estimate  other  rotational  parameters  (e.g.,  the  angular  momentum  vector  or  the  object’s 
moment-of-inertia  ratios)  that  also  generally  require  a  simultaneous  shape  analysis  [19]. 
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4.2  Synodic/Sidereal  Spin  State  Analysis 

Satellite  spin  states  can  be  further  characterized  by  first  detecting  and  then  analyzing  synodic/sidereal  variations,  as 
described  in  detail  in  reference  [2],  Synodic  periods  are  the  apparent  period  measured  from  the  light-curves,  and  can 
be  influenced  by  the  relative  motion  of  the  satellite,  observer,  and  illumination  source.  For  instance,  if  a  spinning 
body  moves  quickly  with  respect  to  an  observer,  the  apparent  rotation  rate  as  seen  from  the  observer’s  perspective 
can  deviate  significantly  from  the  actual  rotation  rate  of  the  body.  (An  extreme  example  in  this  respect  is  the  Moon, 
which  always  presents  the  same  face  to  the  Earth  showing  a  synodic  rotation  rate  of  zero  even  though  it  actually  has 
a  sidereal  rotation  rate  of  1  revolution  every  28  days.)  Although  developed  for  asteroid  characterization,  the 
synodic/sidereal  method  applies  equally  well  to  artificial  satellite  observations  because  it  is  independent  of  the  shape 
and  reflectance  characteristics  of  the  observed  object.  In  fact,  the  sliding  window  PFS  periodicity  analysis  described 
in  section  3.4  was  designed  specifically  to  detect  synodic  periodicity  variations  [2], 

The  synodic/sidereal  method  uses  a  set  of  measured  synodic  frequencies  to  estimate  a  rotating  object’s  sidereal 
angular  spin  rate  and  spin  axis  orientation.  This  method  relies  on  two  fundamental  assumptions: 
1)  that  the  observed  body  is  in  a  stable  spin  rotation  state,  and  2)  that  each  full-cycle  synodic  period  measured  by  the 
observer  corresponds  to  the  time  interval  required  for  the  body  to  rotate  such  that  the  “phase  angle  bisector”  (PAB) 
makes  one  complete  revolution  in  azimuthal  angle  around  the  spin  axis.  The  PAB  is  the  direction  midway  between 
the  body-to-observer  and  body-to-Sun  direction  vectors.  Over  short  time  scales,  corresponding  to  small  changes  in 
observation  and  solar  illumination  angles,  the  single-cycle  reflected-light  signature  of  a  spinning  body  is  controlled 
to  a  large  extent  (but  not  entirely)  by  the  position  of  the  PAB  in  the  body’s  reference  frame.  This  provides  the  basis 
for  the  second  assumption,  and  makes  the  method  shape-independent  to  first  order. 

Synodic/sidereal  spin  state  estimation  analyses  can  be  performed  using  the  algorithm  described  in  reference  [2].  In 
that  case,  the  results  were  shown  using  a  least-squared  chi-squared  metric  plotted  over  the  unit  sphere  using  a 
Mollweide  projection  (see  Figure  7  of  [2]).  The  current  analysis  uses  a  similar  algorithm  to  estimate  a  set  of  three 
parameters,  { 6j,,  t//0,  co},  which  denote  the  two  Euler  angles  of  the  spin  axis,  and  the  sidereal  angular  spin  rate  co  = 
Ays! At  =  2tt.IT,  where  T  is  the  sidereal  full-cycle  period.  The  two  Euler  angles  convert  directly  to  the  inertial  spin 
axis  RA/DEC  angles,  a0  =  </>q  -  ti/2  and  50  =  Jt/2  -  6th  respectively.  The  current  algorithm  also  provides  an  estimate 
for  the  probability  distribution  function  (PDF)  of  the  spin  axis  orientation.  Ideally,  this  PDF  has  a  single  sharp 
maximum,  indicating  that  the  estimation  process  has  converged  to  a  single  a  single,  well-determined  spin  axis.  This 
is  indeed  the  case  when  the  current  method  is  applied  to  the  IMAGE  satellite  data  from  [2],  which  contain  strong, 
well-detected  synodic  periodicity  variations.  However,  in  cases  where  synodic  variations  are  less  strongly  detected, 
the  PDF  can  have  multiple  local  maxima  and/or  be  widely  distributed,  indicating  either  multiple  feasible  solutions 
and/or  large  extended  regions  where  feasible  solutions  exist.  In  these  cases,  improving  the  spin  state  estimate  entails 
acquiring  light -curve  data  that  show  stronger  and/or  better-detected  synodic  variations. 

4.3  Synodic/Sidereal  Spin  State  Analysis  for  Delta-IV  Stage  2  Rocket-  Body  SCN  38834 

The  lower  panel  of  Figure  2  shows  that,  in  this  SCN  38834  data  set,  no  synodic  variations  were  detectable  using  the 
sliding  window  method.  Notably,  this  does  not  mean  that  none  exist,  just  that  none  were  detectable  at  a  sufficient 
level  of  confidence.  This  really  means  that  a  flat  line  fits  the  period  measurements  along  with  their  error  bars,  which 
is  clearly  the  case  shown  in  the  in  the  middle  panel  plot  of  Figure  2.  Synodic  variations  at  a  level  smaller  than  the 
1 -sigma  error  bars  could  (and  probably  do)  exist,  but  cannot  be  detected  in  the  given  data  set.  Notably,  these  error 
bars  can  be  reduced  by  temporal  measurements  with  higher  SNR,  faster  cadence,  and/or  longer  light-curve  duration, 
which  points  to  the  most  natural  means  of  increasing  probabilities  of  synodic  detections  for  this  object.  However, 
given  that  this  this  data  set  contains  no  detectable  synodic  variations,  then  no  synodic/sidereal  spin  state  analysis  can 
be  performed,  and  the  spin  axis  cannot  be  constrained  at  all. 

4.4  Synodic/Sidereal  Spin  State  Analysis  for  Delta-IV  Stage  2  Rocket-  Body  SCN  29156 

The  lower  panel  of  Figure  3  shows  an  example  where  weak  synodic  variations  have  been  detected  for  SCN  29156, 
which,  in  turn,  allows  a  synodic/sidereal  spin  state  estimation  analysis  to  be  performed.  Figure  6  shows  the  resulting 
PDF  of  the  spin  axis  orientation  using  a  Mollweide  projection.  Notably,  in  this  case  the  PDF  is  non-zero  over  a 
large  region  of  the  unit  sphere  has  four  local  maxima,  indicating  multiple  feasible  solutions  —  a  direct  result  of  the 
weakly  detected  synodic  variations. 
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Figure  6.  The  probability  distribution  function  of  the  spin  axis  resulting  from  synodic/sidereal  analysis  of  the 
synodic  periodicities  for  SCN  29156  plotted  in  the  lower  panel  of  Figure  3.  The  widely  distributed,  multi- 
peaked  PDF  indicates  that  the  estimation  process  yielded  multiple  feasible  spin-state  solutions.  In  this  case, 
higher-quality  data  must  be  gathered  to  converge  on  an  isolated,  unique  solution. 


Figure  7.  The  probability  distribution  functions  of  the  spin  axis  resulting  from  two  different  estimation 
analyses  for  the  D4S2  rocket-body  SCN  36586,  both  based  on  the  light-curve  data  shown  in  Figure  4.  The  left 
panel  shows  the  PDF  from  the  shape-independent,  synodic/sidereal  method,  and  the  right  panel  from  the 
shape-dependent,  simple-cylinder-approximation  rotation  state  analysis.  Note  that  both  PDFs  are  plotted  on 
the  same  scale,  and  the  grid  lines  indicate  RA  and  DEC  values  in  the  horizontal  and  vertical  directions 
respectively.  These  sharply-peaked  PDFs  indicate  both  methods  converge  on  well -determined  solutions, 
although  with  somewhat  different  overall  shapes  and  peak  RA/DEC  values. 

4.5  Synodic/Sidereal  Spin  State  Analysis  for  Delta-IV  Stage  2  Rocket-  Body  SCN  36586 

The  lower  panel  of  Figure  4  shows  an  example  where  stronger  synodic  variations  have  been  detected.  Specifically, 
for  that  SCN  36586  2013-12-24  light-curve,  the  sliding -window  method  detects  relatively  strong  periodicity 
variations  in  the  hour  or  so  centered  at  ~  11:30  UT.  Subjecting  these  variations  to  a  synodic/sidereal  analysis  yields 
the  spin  axis  PDF  shown  on  the  left  panel  of  Figure  7.  Notably,  this  PDF  has  a  relatively  sharp  peak,  and  the  figure 
shows  a  Mollweide  projection  of  just  a  portion  of  the  unit  sphere  centered  on  the  peak.  The  best-estimate  spin  state 
in  this  case  has  a  sidereial  period  of  T  =  1 18.2689  ±  0.0051  s,  a  spin  axis  orientation  given  by  J2000  right-ascension 
and  declination  a0  =  99.4°  and  §0  =  27.8°,  with  a  1-sigma  estimation  uncertainty  of  about  1.1°. 


Cleared  for  Public  Release  (Release  #  377ABW-2014-0646,  OPS-14-6494) 


Cleared  for  Public  Release  (Release  #  377ABW-2014-0646,  OPS-14-6494) 


5  SHAPE-DEPENDENT  ROTATION-STATE  ANALYSIS 

The  periodicity-based  methods  described  in  the  previous  section  provide  rotation  state  information  in  a  shape - 
independent  fashion,  i.e.,  without  any  explicit  consideration  of  the  object’s  shape  or  reflectance  characteristics. 
Alternatively,  shape-dependent  rotation  state  estimation  processes  can  be  formulated.  These  come  in  three  basic 
categories:  1)  those  in  which  all  shape/reflectance  information  is  estimated  from  the  data  along  with  the  rotation 
state  information),  2),  those  where  some  but  not  all  of  the  shape/reflectance  information  is  estimated,  or  3)  those 
where  no  shape/reflectance  information  is  estimated  (because  it  is  all  known  or  assumed  a  priori).  Most  shape- 
dependent  asteroid  methods  fall  into  the  second,  hybrid  category  (see  [17-19]  and  references  therein),  because  these 
methods  usually  assume  some  a  priori  knowledge  of  the  asteroid  surface  bi-directional  reflectance  distribution 
function  (BRDF),  and  then  estimate  the  asteroid’s  shape  via  a  best-fit  parameterized  model.  This  hybrid  approach 
can  be  used  because  asteroid  BRDFs  are  relatively  well  studied,  but  their  shapes  are  generally  unknown  [17], 
Asteroid  shape  models  range  considerably  in  complexity:  from  simple  (tri-axial  ellipsoids),  moderate  (spherical- 
harmonic  surface  expansion,  multi-vertex  convex-hulls),  up  to  relatively  complex  (multi-vertex  surfaces  with 
concavities). 

5.1  Rotation  State  Estimation  for  Satellites  of  Known  Shape  and  Reflectance 

For  man-made  satellites  with  known  shapes  and  materials,  one  may  be  tempted  to  evaluate  rotation  states  using  the 
third  method  discussed  above,  in  which  no  shape/reflectance  information  about  the  object  is  estimated  in  the 
analysis  process.  This  estimation  process  requires  having  a  complete  model  of  the  shape  and  reflectance 
characteristics  of  a  satellite,  and  then  using  a  forward-modeling  process  such  as  that  described  in  [20]  in  order  to 
estimate  only  the  object’s  attitude  parameters,  an  approach  used  in  [21]  for  stabilized  satellites  and  [22]  for  satellites 
that  have  potentially  maneuvering  attitudes.  For  objects  in  torque-free  multi-axis  rotation  states,  this  kind  of 
analysis  would  entail  estimating  the  following  set  of  eight  parameters  [18,  19]:  {$ },  £j>,  y/0,  L,  i\,  i2}.  Here  ( ,  Oo, 
i/At)  denote  the  body’s  orientation  Euler  angles  for  a  given  initial  time,  f0,  L  is  the  3-component  angular  momentum 
vector,  and  (q,  i2)  are  the  two  smallest  principal-axis-frame  moments  of  inertia,  normalized  to  the  largest  moment  of 
inertial  (i.e.,  q  =  /,//3,  and  /j  =  /2//3).  In  this  case  the  rotation  state  estimation  process  entails  solving  an  eight¬ 
dimensional  optimization  problem  [19].  However,  this  dimensionality  can  be  reduced  by  one  for  bodies  with  axi- 
symmetric  mass  distributions  because  ;j  =  i2,  a  good  approximation  for  many  asteroids  [17]  and  man-made  rocket- 
bodies  [23],  Furthermore,  for  objects  in  torque-free  principal-axis  spin  states,  the  equations  of  motion  no  longer 
depend  on  (;j,  i2 ),  and  the  set  of  estimated  parameters  can  be  reduced  to  four,  { (f>0,  do,  y/o,  co},  where,  as  before,  the 
first  two  Euler  angles  provide  the  spin  axis  orientation  and  co  denotes  the  angular  spin  rate. 

The  rotation  state  parameters  outlined  above  strictly  apply  to  objects  in  torque-free  environments.  However, 
satellites  in  Earth  orbit  are  perturbed  by  environmental  torques  [e.g.,  23].  In  this  case,  the  parameters  outlined  above 
would  still  apply,  but  must  be  estimated  from  light-curve  data  collected  over  a  time-scale  sufficiently  short  so  that 
perturbations  from  on-orbit  torques  can  be  safely  neglected. 

5.2  Rotation  State  Estimation  for  Objects  of  Unknown  Shape  and  Unknown  Reflectance  Characteristics 

Reference  [24]  outlines  an  attitude  evaluation  method  requiring  no  a  priori  shape  or  reflectance  information  (also 
see  [25]).  In  this  case,  the  surface  of  the  man-made  satellite  being  observed  is  modeled  as  an  ensemble  of  flat  facets, 
and  the  size  and  reflective  properties  of  each  facet  estimated  from  the  light-curve  data  along  with  the  rotation  state 
parameters.  The  sizes  of  the  facets  can  be  modeled  in  terms  of  albedo-area  products  [24]  or  as  projected  areas  [25]. 
Different  surface  compositions  on  the  satellite’s  surfaces  can  be  accounted  for  by  estimating  different  specular  vs. 
diffuse  components  for  individual  facets  [24],  adjusted  to  match  shiny  vs.  dull  materials.  This  approach  has  the 
advantage  of  being  applicable  to  completely  unknown  satellites,  but  increases  the  dimensionality  of  the  optimization 
analysis  considerably,  sometimes  to  a  computationally-prohibitive  extent. 

5.3  Rotation  State  Estimation  for  Objects  of  Known  Shape  but  Unknown  Reflectance  Characteristics 

Often  satellites  have  known  shapes,  but  the  materials  covering  their  surfaces  are  unknown  or  have  uncharacterized 
reflectance  characteristics.  In  these  cases,  a  hybrid  approach  is  possible  where  only  the  surface  BRDF  parameters 
are  estimated  along  with  the  rotation  state  parameters.  These  BRDF  parameters  could  comprise  laboratory - 
measured  material  reflectance  functions,  in  which  case  the  process  would  entail  estimating  which  materials  on 
which  satellite  surfaces  (or  components)  provide  the  best  match  to  the  observed  data,  which  might  be  formulated  in 
a  manner  similar  to  that  described  in  [26],  or  as  a  linear-programming  or  discrete  goal-maximization  optimization 
analysis.  Alternatively,  the  parameters  could  comprise  a  set  of  analytical  BRDF  expansion  functions  such  as  those 
described  in  [27],  which  would  naturally  adjust  to  match  surface  material  shiny  vs.  dull  reflectance  properties,  which 
is  the  approach  used  in  this  analysis.  The  known  shapes  could  comprise  wireframe  models  or  solid  models  for 
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complex  satellites,  or  more  basic  shapes  to  approximate  simpler  objects  (e.g.,  cubes  or  for  CubeSats,  cylinders  for 
rocket-bodies,  etc.). 

5.4  Rotation  State  Estimation  for  Nearly  Cylindrical  Objects 

Some  satellites  can  be  approximated  reasonably  well  using  simple  cylindrical  shapes  to  model  reflected 
brightnesses.  These  include  Boeing  376  GEO  satellites  —  dual-spinning  cylindrical  satellites  with  outer  round 
surface  covered  in  shiny  solar  cells  —  which  can  produce  bright,  short-duration  flashes  when  they  glint  sunlight  [1, 
4],  Because  the  timing  of  the  bright  flashes  can  be  measured  very  accurately,  synodic  periodicity  variations  for 
these  objects  can  be  tracked  much  better  than  for  less  shiny  objects. 

Upper-stage  rocket-bodies  can  also  be  approximated  as  cylinders,  although  perhaps  less  accurately  than  Boeing  376 
satellites.  Figure  1  shows  that  the  D4S2  rocket-bodies  do  have  a  roughly  cylindrical  shape,  and  much  of  their 
projected  area  comprises  a  nearly  perfectly  cylindrical  white-colored  surface.  Our  observations  indicate  that  the 
specular  component  of  reflection  for  D4S2s  is  very  low,  which  means  that  they  generally  do  not  produce  bright, 
short-duration  flashes.  This  is  demonstrated  in  the  bottom  panel  of  Figure  5,  which  shows  an  observed  D4S2 
rotational  light-curve  characterized  by  the  two-peaks-per-cycle  brightness  pattern  expected  for  a  cylinder  in 
principal  axis  rotation  [28],  Notably  these  brightness  peaks  are  not  short,  bright  spikes,  as  is  observed  for  Boeing 
376  satellites  [1,  4],  Because  of  this,  non-specular  rocket-bodies  such  as  D4S2s  do  not  lend  themselves  as  readily  to 
synodic/sidereal  spin  state  analyses.  This  is  especially  true  for  rocket-bodies  in  orbits  that  have  high  perigee 
altitudes  (e.g.,  GPS  D4S2s),  that  move  at  very  low-angular  rates  relative  to  ground-based  observers,  which  also 
makes  synodic  periodicity  variations  more  difficult  to  detect  [2],  In  fact,  multiple  Raven  light-curve  observations  of 
the  D4S2s  that  have  launched  GPS  satellites  have  been  conducted,  and  most  yield  no  detection  of  synodic  variations 
(as  shown  in  the  bottom  panel  of  Figure  2)  or  only  weak  detections  (Figure  3,  bottom  panel).  In  fact,  the  synodic 
variations  shown  in  the  bottom  panel  of  Figure  4  represent  the  strongest  detected  in  about  a  dozen  D4S2  light-curves 
of  comparable  quality.  Our  analysis  indicates  that  those  variations  were  likely  detected  only  because  during  that 
light-curve  the  PAB  vector  made  a  relatively  close  approach  to  the  spin  axis,  an  effect  which  tends  to  create  large 
synodic  variations  [1,2], 

These  difficulties  of  measuring  synodic  variations  for  non-specular  rocket-bodies  prompt  an  investigation  of  a 
shape-dependent  spin-state  estimation  process.  For  this,  brightnesses  are  approximated  using  a  simple  cylinder 
shape  model,  characterized  by  flat  top  and  bottom  surfaces  with  pine  Lambertian  reflectances,  and  a  round  surface 
with  Lambertian  +  specular  reflectances.  Each  surface  is  parameterized  using  albedo-area  (aA)  products  (see  [24, 
25]).  The  top  and  bottom  surfaces  are  parameterized  by  ( aAT ,  aAn).  The  round  surface  is  parameterized  by  a 
Lambertian  albedo-area  product,  aAs,  as  well  as  a  series  of  specular  functions  aAm,  where  the  index  m  spans  a  range 
of  specular  functions  such  as  a  series  of  Phong  or  Cook-Torrance  BRDFs  (see  [24]  for  more  details  and  examples). 
These  aA  products  must  be  estimated  along  with  the  rotation  state  parameters  using  a  process  that  seeks  the  best  fit 
to  the  observed  light-curve  data. 

5.5  Spin  State  Estimation  for  Delta-IV  Stage  2  Rocket-bodies 

For  this  study  we’ve  developed  and  implemented  a  shape-dependent  algorithm  to  estimate  spin  states  using  the 
simple-cylinder  formulation  outlined  above.  Our  analysis  of  multiple  D4S2  data  sets  indicates  that  adding  specular 
components  to  the  shape  model  does  not  improve  the  quality  of  the  fits  to  the  light-curves  in  a  statistically  justifiable 
way.  So  for  D4S2s,  seven  parameters  must  be  estimated  for  the  shape-dependent  spin  state  analysis:  { </;0,  6j,,  t//0,  co, 
ciAj,  «4|i,  a  A  s } .  While  still  numerically  tractable,  this  is  significantly  larger  than  the  three  required  for  the  shape- 
independent  synodic/sidereal  analysis,  { (j>0 ,  90,  co],  as  described  in  section  4.2.  Space  limitations  prevent  a  detailed 
description  of  our  newly-implemented  shape -dependent  formulation  in  this  paper,  but  a  few  salient  features  based  on 
analyses  of  simulated  rocket-body  data  sets  warrant  comment: 

1.  The  seven-dimensional  parameter  space  for  the  estimation  process  is  characterized  by  highly  non-linear 
objective  functions  with  multiple  local  extrema  (some  cyclical  in  nature),  so  confidently  finding  a  global 
solution  to  the  optimization  problem  can  be  exceedingly  computationally  intensive. 

2.  The  search  range  for  the  spin  rate  co  should  be  limited  to  a  realistic  span  provided  by  a  previously-calculated 
periodicity  analysis,  allowing  for  some  reduction  of  computation. 

3.  Given  a  hypothesized  spin  state  { $>,  y/o,  co],  best  estimates  for  the  aA  products  can  be  found  in  a 
straightforward  manner  [2],  also  providing  a  means  to  reduce  computation. 

4.  Employing  data  with  phase  angles  greater  than  about  60°  to  90°  can  lead  to  both  biases  and  inaccuracies  in  the 
estimation  process,  apparently  because  the  simple-cylinder  shape  model  approximation  fails  to  match  both  high- 
and  low-phase-angle  data  simultaneously. 
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5.  Our  current  analytical  formulation  for  uncertainties  and  covariances  tends  to  significantly  underestimate 
deviations  from  known  truth  for  the  simulated  rocket -body  analyses.  This  also  may  be  due  to  the  oversimplified 
shape  model,  but  warrants  further  study.  Numerical  boot-strap  or  jack-knife  uncertainty  estimates  would  likely 
provide  better  agreement,  but  are  computationally  intensive  and  have  not  been  implemented  for  this  study. 

6.  Although  the  shape -dependent  formulation  does  not  require  that  synodic  variations  be  present  in  the  input  data, 
it  turns  out  that  when  none  exist,  the  analysis  generally  provides  two  equal-quality,  degenerate  spin  state 
solutions  with  oppositely-directed  spin  axis  orientations.  In  other  words,  the  shape-dependent  method  has 
difficultly  discriminating  between  clockwise  and  counter-clockwise  spins.  This  result  is  not  surprising,  and  can 
be  predicted  from  simple  symmetry  considerations.  However,  given  light-curves  that  do  include  strong  synodic 
variations,  a  properly  formulated  shape-dependent  estimation  analysis  should  be  able  to  converge  to  a  unique 
spin  axis,  just  like  the  synodic/sidereal  method  does. 

The  right  panel  of  Figure  7  shows  the  spin  axis  PDF  calculated  by  applying  the  shape -dependent  estimation  process 
describe  above  to  the  D4S2  SCN  36586  2013-12-24  light-curve  shown  in  the  top  panel  of  Figure  4,  but  restricted  to 
phase  angles  <  70°  to  limit  bias  and  inaccuracies.  Note  that  this  shape-dependent  PDF  is  plotted  on  the  exact  same 
scale  as  the  shape-independent  PDF  shown  in  the  left  panel,  and  is  based  on  the  same  data  set.  The  shape-dependent 
analysis  also  yields  a  PDF  with  a  relatively  sharp  peak,  but  a  somewhat  different  best-estimate  spin  state:  sidereal 
period  T  =  118.2940  ±  0.0036  s,  and  spin  axis  J2000  RA/DEC  values  a0  =  103.8°  and  50  =  24.0°,  with  a 
analytically-estimated  1 -sigma  estimation  uncertainty  of  about  0.2°.  These  parameters  agree  reasonably  well  with 
the  synodic/sidereal  values  (reported  in  section  4.5):  the  sidereal  periods  only  differ  by  ~  0.025  s  and  the  spin  axis 
unit  vectors  by  ~  5.5°.  However,  these  differences  significantly  exceed  the  estimated  1-sigma  ranges  for  both 
methods,  which  we  believe  to  be  due  (at  least  in  part)  to  systematic  under-estimation  of  uncertainties  that  the  shape- 
dependent  analysis  provides  (at  least  in  its  current  formulation).  The  process  also  yields  aA  estimates  of  <  1.5  m2 
for  the  two  cylinder  ends  and  ~  42  m2  for  the  round  side,  consistent  with  the  known  dimensions  and  materials  for 
D4S2  rocket-bodies. 

It  is  also  worth  noting  that  calculating  the  shape-dependent  PDF  shown  on  the  right  panel  of  Figure  7  required 
approximately  a  factor  of  100  more  computation  time  than  the  shape -independent  synodic/sidereal  PDF  on  the  left 
panel.  This  difference  arises  because  the  shape -dependent  analysis  has  significantly  higher  dimensionality  and  non¬ 
linearity  than  the  shape-independent  approach.  It  is  likely  that  for  more  complicated  satellites,  this  computational 
difference  will  be  even  more  dramatic  —  because  those  objects  will  require  higher  dimensionality  shape  models 
than  the  simple  cylinder  employed  here  —  probably  requiring  the  use  of  parallel,  multi-node  computing. 

6  CONCLUSIONS  AND  FUTURE  WORK 

Our  study  of  D4S2  rocket-bodies  satellites  indicate  that  synodic  variations  can  be  difficult  to  detect  from  the  ground 
for  these  non-specular,  nearly-cylindrical  objects.  Multiple  Raven  D4S2  light-curve  observations  were  conducted 
for  this  study,  with  most  yielding  none  or  weak  synodic  detections.  Only  one  light-curve  provided  relatively  strong 
synodic  variations,  for  SCN  36586  as  shown  in  Figure  4.  Spin  state  analyses  of  this  light-curve  using  the  shape- 
independent  synodic/sidereal  approach  as  well  as  a  shape-dependent  simple-cylinder  approach  yield  roughly  the 
same  results  for  the  object’s  spin  state  parameters  —  although  significantly  outside  the  estimated  1-sigma  limits  of 
one  another,  which  warrants  more  investigation.  So  our  analysis  indicates  that  a  simple-cylinder  shape  model  can,  in 
principle,  suffice  to  characterize  rotation  states  for  upper-stage  rocket  bodies  observed  at  low  phase  angles,  although 
more  work  needs  to  be  accomplished. 

Overall,  the  difficulty  in  detecting  synodic  variations  from  the  ground  for  objects  in  high-altitude  orbits  compels  an 
investigation  of  shape-dependent  spin  state  analysis  methods,  which  do  not  require  detection  of  synodic  variations. 
Unfortunately,  shape-dependent  approaches  seem  be  relatively  computationally  intensive,  and  often  yield  two 
degenerate  solutions  with  oppositely-directed  spin  axes.  For  these  reasons,  the  shape-independent  synodic/sidereal 
spin-state  estimation  approach  should  be  favored  whenever  feasible,  which  sometimes  entails  acquiring  light-curve 
data  that  show  stronger  and/or  better-detected  synodic  variations.  In  this  regard,  it  is  worth  noting  that  even  for 
deep-space  objects,  one  can  significantly  improve  the  detectability  of  synodic  variations  even  for  non-specular 
objects  by  using  a  space-based  sensor  [29], 
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